FUNCTION IReqw, lambda, flux, linecent, linewidth, cont1, cont1width, cont2, cont2width, PLOT = PLOT

;interpolate the continuum over the line
incont1 = WHERE(lambda GT cont1-cont1width/2. AND lambda LT cont1+cont1width/2., n_incont1)
incont2 = WHERE(lambda GT cont2-cont2width/2. AND lambda LT cont2+cont2width/2., n_incont2)
inline = WHERE(lambda GT linecent-linewidth/2. AND lambda LT linecent+linewidth/2., n_inline)

IF n_incont1 GT 1 AND n_incont2 GT 1 AND n_inline GT 1 THEN BEGIN

    ;find the median flux in the continuum regions
    mediancont1 = MEDIAN(flux[incont1])
    mediancont2 = MEDIAN(flux[incont2])
    
    ;interpolate between continuum
    ;regions through the line
    continuum = INTERPOL([mediancont1,mediancont2],[cont1,cont2],lambda[inline])

    ;find the flux absorbed in the line
    absorbed_flux = TSUM(lambda[inline], continuum - flux[inline])

    ;find the mean flux in the line
    line_flux = TSUM(lambda[inline], flux[inline])/(lambda[inline[n_inline-1]]-lambda[inline[0]])

    ;divide by the continuum flux density in the center of the line
    distance_to_linecenter = linecent-lambda[inline]
    closest = MIN(ABS(distance_to_linecenter), closest_index)
    eqw = absorbed_flux*1./continuum[closest_index]
    flux_ratio = line_flux*1./continuum[closest_index]

    ;make the index a magnitude
    index = -2.5*ALOG10(flux_ratio)

    IF KEYWORD_SET(PLOT) THEN BEGIN
        plotstart = cont1-4*cont1width
        plotend = cont2+4*cont2width

        onplot = WHERE(lambda GT plotstart AND lambda LT plotend)
        hiplot = MAX(flux[onplot])
        loplot = MIN(flux[onplot])

        PLOT, lambda, flux, XRANGE = [plotstart,plotend],YRANGE = [0.8*loplot,1.2*hiplot], PSYM = 10, /XSTY, /YSTY
        OPLOT, lambda[inline], continuum, LINESTYLE = 2
        OPLOT, [linecent+linewidth/2.,linecent+linewidth/2.], [-1000.,1000.], LINESTYLE = 3
        OPLOT, [linecent-linewidth/2.,linecent-linewidth/2.], [-1000.,1000.], LINESTYLE = 3
        OPLOT, [cont1-cont1width/2.,cont1-cont1width/2.], [-1000.,1000.], LINESTYLE = 4
        OPLOT, [cont1+cont1width/2.,cont1+cont1width/2.], [-1000.,1000.], LINESTYLE = 4
        OPLOT, [cont2-cont2width/2.,cont2-cont2width/2.], [-1000.,1000.], LINESTYLE = 4
        OPLOT, [cont2+cont2width/2.,cont2+cont2width/2.], [-1000.,1000.], LINESTYLE = 4
        OPLOT, [lambda[inline[closest_index]],lambda[inline[closest_index]]], [continuum[closest_index], continuum[closest_index]], PSYM = 4
        
        wait = GET_KBRD(1)

    ENDIF

ENDIF ELSE BEGIN
    eqw = -99.
    index = -99.
ENDELSE

RETURN, eqw*10000. ;index ;

END
;-----------------------------------------------
FUNCTION IRratio, lambda, flux, numcent, numwidth, denomcent, denomwidth, PLOT = PLOT

innum = WHERE(lambda GT numcent-numwidth/2. AND lambda LT numcent+numwidth/2., n_num)
indenom = WHERE(lambda GT denomcent-denomwidth/2. AND lambda LT denomcent+denomwidth/2., n_denom)

IF n_num GT 1 AND n_denom GT 1 THEN BEGIN

    ratio = MEAN(flux[innum]/flux[indenom])

    IF KEYWORD_SET(PLOT) THEN BEGIN
        plotrange = MAX([lambda[innum],lambda[indenom]]) - MIN([lambda[innum],lambda[indenom]])
        plotstart = MIN([lambda[innum],lambda[indenom]]) - plotrange/2.
        plotend = MAX([lambda[innum],lambda[indenom]]) + plotrange/2.

        onplot = WHERE(lambda GT plotstart AND lambda LT plotend)
        hiplot = MAX(flux[onplot])
        loplot = MIN(flux[onplot])

        PLOT, lambda, flux, XRANGE = [plotstart,plotend],YRANGE = [0.8*loplot,1.2*hiplot], PSYM = 10, /XSTY, /YSTY
        OPLOT, [numcent-numwidth/2.,numcent-numwidth/2.], [-1000.,1000.], LINESTYLE = 4
        OPLOT, [numcent+numwidth/2.,numcent+numwidth/2.], [-1000.,1000.], LINESTYLE = 4
        OPLOT, [denomcent-denomwidth/2.,denomcent-denomwidth/2.], [-1000.,1000.], LINESTYLE = 4
        OPLOT, [denomcent+denomwidth/2.,denomcent+denomwidth/2.], [-1000.,1000.], LINESTYLE = 4
        
        wait = GET_KBRD(1)

    ENDIF

ENDIF ELSE BEGIN

    ratio = -99.

ENDELSE

RETURN, ratio

END

;-----------------------------------------------
FUNCTION measureIRindices, lambda, flux
 
;define an array where I can store the spectral indices
indices = FLTARR(30)

;start measuring indices

;measure Mg I @ 1.48 microns in the H band
indices[0] = IREQW(lambda, flux, 1.4892, 0.006, 1.482, 0.003, 1.494, 0.003)

;measure Mg I @ 1.50 microns in the H band
indices[1] = IREQW(lambda, flux, 1.504, 0.004, 1.498, 0.0045, 1.5095, 0.0045)

;measure K I @ 1.52 microns in the H band
indices[2] = IREQW(lambda, flux, 1.5172, 0.004, 1.5105, 0.004, 1.523, 0.004)

;measure Mg I @ 1.576 
indices[3] = IREQW(lambda, flux, 1.576, 0.004, 1.57, 0.004, 1.58, 0.004)

;measure a simple ratio using two equivalent widths with the same
;continuum regions...
indices[4] = -IREQW(lambda, flux, 1.575, 0.002, 1.57, 0.004, 1.5815, 0.003)+IREQW(lambda, flux, 1.5775, 0.003, 1.57, 0.004, 1.5815, 0.003)
;indices[4] = IRRATIO(lambda, flux, 1.575, 0.002, 1.5775, 0.003, /plot)

;measure Si I @ 1.59 microns in the H band
indices[5] = IREQW(lambda, flux, 1.59, 0.005, 1.586, 0.003, 1.594, 0.003)

;measure CO @ 1.62 microns in the H band
indices[6] = IREQW(lambda, flux, 1.61975, 0.0045, 1.615, 0.003, 1.628, 0.003)

;measure CO @ 1.66 (inspired by Meyer et al. 1998)
indices[7] = IREQW(lambda, flux, 1.6625, 0.004, 1.6565, 0.003, 1.67, 0.003)

;measure Al @ 1.67 microns in the H band
indices[8] = IREQW(lambda, flux, 1.674, 0.006, 1.659, 0.003, 1.678, 0.002)

;measure HI @ 1.68 microns in the H band
indices[9] = IREQW(lambda, flux, 1.681, 0.004, 1.6775, 0.003, 1.6855, 0.003)

;measure OH in H band (inspired by Meyer et al. 1998)
indices[10] = IREQW(lambda, flux, 1.686, 0.0065, 1.6785, 0.003, 1.695, 0.003)

;Measure bump next to 1.71 Mg line
indices[11] = IREQW(lambda, flux, 1.7075, 0.003, 1.704, 0.003, 1.7145, 0.003)

;Measure Mg I @ 1.71 microns in the H band
indices[12] = IREQW(lambda, flux, 1.7115, 0.003, 1.704, 0.003, 1.7145, 0.003)

;measure a continuum slope in the H band
indices[13] = IRRATIO(lambda, flux, 1.605, 0.02, 1.69, 0.02)

;measure water absorption in the H band
indices[14] = IRRATIO(lambda, flux, 1.69, 0.02, 1.77, 0.02)

;Measure Ca I @ 1.982 microns in the H band
indices[15] = IREQW(lambda, flux, 1.982, 0.013, 1.9676, 0.005, 1.99775, 0.005)

;measure Mg I @ 2.1066 + Al I 2.1099 microns in the K band 
indices[16] = IREQW( lambda, flux, 2.1082, 0.006, 2.102, 0.004, 2.113, 0.004)

;measure Mg I @ 2.1066 microns in the K band 
indices[17] = IREQW( lambda, flux, 2.1065, 0.003, 2.102, 0.004, 2.113, 0.004)

;measure Al I @ 2.1099 microns in the K band 
indices[18] = IREQW( lambda, flux, 2.1097, 0.003, 2.102, 0.004, 2.113, 0.004)

;measure Al I @ 2.1168 microns in the K band
indices[19] = IREQW( lambda, flux, 2.1168, 0.003, 2.1145, 0.002, 2.1195, 0.002)

;measure Brackett Gamma @ 2.165 microns in the K band (using
;bandpasses from Kleinmann & Hall 86, but with a modified red continuum)
indices[20] = IREQW( lambda, flux, 2.16625, 0.0047, 2.155, 0.0044, 2.175, 0.0042)

;measure Na I doublet @ 2.21 microns in the K band (using bandpasses
;from Ali et al. 95)
;indices[21] = IREQW( lambda, flux, 2.2075, 0.007, 2.197, 0.003 , 2.2185, 0.003, /plot)
indices[21] = IREQW( lambda, flux, 2.2075, 0.007, 2.195, 0.004, 2.217, 0.006)

;measure asymmetry in NaI profile for giant/dwarf separation
indices[22] = IREQW( lambda, flux, 2.206, 0.0035, 2.195, 0.004, 2.217, 0.006)
indices[23] = IREQW( lambda, flux, 2.2095, 0.0035, 2.195, 0.004, 2.217, 0.006);/IREQW( lambda, flux, 2.206, 0.0035, 2.195, 0.004, 2.217, 0.006)

;measure Ca I @ 2.26 microns in the K band (using bandpasses
;from Ali et al. 95)
indices[24] = IREQW( lambda, flux, 2.264, 0.007, 2.258, 0.0045, 2.27, 0.0045)

;measure Mg I @ 2.28 microns in the K band (using bandpasses
;from Ivanov et al. 04)
indices[25] = IREQW( lambda, flux, 2.2806, 0.0052, 2.275, 0.006, 2.2857, 0.005)

;measure CO @ 2.3 microns in the K band (using bandpasses 
;from Ivanov et al. 06, plus a long wavelength anchor I invented)
indices[26] = IREQW(lambda, flux, 2.30375, 0.0225, 2.288, 0.007, 2.3185, 0.004)

;measure a continuum slope in the K band
indices[27] = IRRATIO(lambda, flux, 2.19, 0.02, 2.28, 0.02)

;measure water absorption in the K band
indices[28] = IRRATIO(lambda, flux, 2.28, 0.02, 2.37, 0.02)

;measure co feature at 2.345 microns
indices[29] = IREQW(lambda, flux, 2.3455, 0.003, 2.3425, 0.003, 2.349, 0.003)

RETURN, indices

END
